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We present computer simulations of the response of a flexoelectric blue phase network, either in bulk or under confinement, to 
an applied field. We find a transition in the bulk between the blue phase I disclination network and a parallel array of disclinations 
along the direction of the applied field. Upon switching off the field, the system is unable to reconstruct the original blue phase 
but gets stuck in a metastable phase. Blue phase II is comparatively much less affected by the field. In confined samples, 
^ the anchoring at the walls and the geometry of the device lead to the stabilisation of further structures, including field-aligned 
V/~) disclination loops, splayed nematic patterns, and yet more metastable states. Our results are relevant to the understanding of the 
CN switching dynamics for a class of new, "superstable", blue phases which are composed of bimesogenic liquid crystals, as these 
materials combine anomalously large flexoelectric coefficients, and low or near-zero dielectric anisotropy. 
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1 Introduction 

The blue phases (BPs) of cholesteric liquid crystals provide 
a remarkable example of soft materials with a rare combina- 
tion of highly nontrivial geometric and topological nature, and 
promising technological potential. BPs are three-dimensional 
networks of disclination lines, which appear due to the spon- 
taneous tendency of a cholesteric liquid crystal to twist along 
i more than one direction. The associated double twist re- 
gions lead to topological frustration, as they cannot be used 
to smoothly tile the three-dimensional space, so that defects, 
or disclination lines, are unavoidable. In two of the experi- 
mentally observed blue phases, BP I and BP II, to which we 
restrict the present study, such disclinations form ordered lat- 
tices, with symmetries characterised by point groups normally 
with associated crystals (BP III is instead likely amorphous). 
The spatial period of the disclination lattices is in the submi- 
crometer range, which endows them with interesting optical 
properties - for instance, a swelling or shrinking of the unit 
cell would normally be linked to a change in the colour in the 
sample. 

A leap forward in blue phase based technology has been 
made possible by the stabilisation of these disclination net- 
works over a range of over 60 K, including room tempera- 
ture^ - -. This was a dramatic advance with respect to the pre- 
vious framework, in which BPs were typically only stable over 
a fraction of a degree. As a consequence, it has now become 
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possible to think of blue phases as a possible basis for de- 
vice or display technology, and this has culminated in the re- 
cent fabrication of a blue-phase display device with very fast 
switching and response times 6 . Such switching times are a 
natural consequence of the small scale of the unit cells which 
need to rearrange when a field is applied. Another potential 
technological advantage of blue phases, which has not been 
exploited much in practice so far, is that they are associated 
with a complex and rich free energy landscape, with several 
competing equilibria which are ideal to build bistable or mul- 
tistable devices. Such devices can retain memory of two or 
more states even after the field used to reach these has been 
switched off (this is a highly sought after property in devices 
as it leads to huge savings in energy consumption) 7 . 

While the stable BPs of Ref . 1 involved dispersing polymers 
inside BPs, those of Ref.-, which provide the main motivation 
for our current work, were a one-component fluid of "bimeso- 
genic" Q liquid crystals, which are associated with very large 
flexoelectric coefficients 8 . Flexoelectricity in liquid crystals 
is a linear coupling between applied electric field and induced 
distortion. This arises, for instance, because the microscopic 
shape of the molecules (which can be e.g. pear-shaped) al- 
lows a coupling between splay (or bend) and polarisation. 
The treatment of Ref. 9 shows that in flexoelectric materials 
this coupling leads to a renormalisation of the elastic con- 
stants, which effectively decrease and facilitates the forma- 
tion of disclination lines hence of BPs. Conversely, in ma- 
terials with large flexoelectric coefficients, elastic distortions 
readily can be created by an external electric field - the flex- 
oelectric coupling is a linear coupling between electric field 
and order parameter variations, which is allowed in the the- 

* To form a bimesogenic liquid crystal, one needs a compound which is made 
up by two polar mesogenic units linked by a flexible spacer 8 . 



ory due to the nonsymmetric nature of the molecule. For low 
enough voltage, this contribution will dominate over the fa- 
miliar, quadratic, dielectric coupling between order parameter 
and electric field. 

Because any-device oriented application of blue phases is 
likely to use the stabilised versions of these materials, it is very 
important to understand the switching behaviour of BPs in re- 
sponse to a flexoelectric, rather than dieletric, coupling to an 
electric field, as this is likely to be the dominant contribution 
at least when preparing BP samples as in Ref . 3 . Therefore, 
in this work our goal is to elucidate the effect of flexoelectric- 
ity in the switching of BP domains, or confined BPs, which 
can now be potentially used to build real devices. To this end, 
we use large scale three-dimensional lattice Boltzmann simu- 
lations, which have already allowed in the past a careful study 
of the thermodynamic and flow properties of BPs— - —. Im- 
portantly, our approach correctly incorporates hydrodynamic 
interactions which are known to affect the switching pathway 
of liquid crystal devices. Our work provides a direct sys- 
tematic numerical investigation of the switching dynamics of 
cholesteric blue phases (both in bulk and under confinement) 
in response to a flexoelectric coupling with the applied electric 
field (we refer to this phenomenon as "flexoelectric switch- 
ing"). 

This paper is structured as follows. In Section II, we review 
the free energy density on which the continuum model we use 
is built, taking care to derive and discuss the extra terms which 
are due to flexoelectricity and which will play a key role in our 
work. In Section III we present our results, first focussing on 
the effect of bulk flexoelectricity on the unit cells of a large BP 
sample, and then analysing the case of confined geometries. 
Finally, we draw our conclusion in Section IV, where we also 
compare the flexoelectric switching studied here with previous 
works studying response of BPs to a dielectric coupling. 

2 Models and methods 

2.1 Free energy and equations of motion 

In order to investigate the switching properties of flexoelec- 
tric blue phases, we consider the formulation of liquid-crystal 
hydrodyamics given by Beris and Edwards 15 , generalized for 
cholesteric liquid crystals. In this formulation the equations 
of motion are written in terms of a tensorial order parameter, 
Q. The latter is related to the local direction of the individ- 
ual molecules, n, by Q a $ = (n a n^ — ^§ a |3)> where the angu- 
lar bracket denotes a coarse-grained average and the Greek 
indices label the Cartesian components of Q. The tensor Q 
is traceless and symmetric and its largest eigenvalue, 2/3q, 
(0 < q < 1), measures the magnitude of local order. The equi- 
librium properties of the system are described by a Landau-de 
Gennes free energy density—"—. This comprises a bulk term 



(summation over repeated indices is implied hereafter), 

h = y (1 - |)f£p " ^CapCpyCya + ^(G^p) 2 , (D 

which describes a first-order transition from the isotropic to 
the chiral (cholesteric) phase, together with an elastic contri- 
bution which for cholesterics is— 

fd = f (9 P e«p) 2 + f (e^8^e 8 p+ — Gap) 2 - (2) 
z z po 

In the formula above K is the elastic constant (here we are 
assuming the one-elastic constant approximation) and po is 
the intrinsic helix pitch of the cholesteric. The tensor £ a ^§ is 
the Levi-Civita antisymmetric third-rank tensor, Ao is a con- 
stant (with units of pressure) and y controls the magnitude of 
order (it plays the role of an effective temperature or concen- 
tration according to whether the cholesteric liquid crystal is 
thermotropic or lyotropic). 

To focus on the effect of flexoelectric coupling with an ex- 
ternal field E a , we assume that the liquid crystal has zero di- 
electric anisotropy and the electric field is uniform throughout 
the sample. Flexoelectric contributions can be divided into a 
bulk and a surface one. To account for the former, we add the 
following term to the free energy density 19 : 

/e = £fiQa$( E ady — £y3 a )<2py, (3) 

where > is the bulk flexoelectric constant. Note that each 
term of the flexoelectric coupling is linear in E a and quadratic 
in Q. This is different from the standard dielectric coupling 
where the field couples linearly to the orientational order pa- 
rameter and quadratically with the electric field. The reduced 
potential is V z = E Z L Z , where E z is the component of the elec- 
tric field along the z-direction and L z the size of the sample in 
the same direction. 

To describe the anchoring of the director field and the flex- 
oelectric interactions at the boundary surfaces, the following 
surface contributions are added to the free energy density— 

f s = ^{Q^-<^? + \[Q*^ 

+e} / (a p e ap )£ a . (4) 

The first term is a quadratic contribution ensuring a soft pin- 
ning of the director field on the boundary surface, along a 
direction no. The parameter Wo controls the strength of the 
anchoring and <2 a y = Soften® — Say/3), where So determines 
the magnitude of surface order. In the strong anchoring regime 
Wo is practically infinite so that the director field is effectively 
fixed at the boundaries. In the weak anchoring regime a term 
proportional to K, describing the elastic distorsion at the sur- 
faces, is also important (see second term in Eq. . Finally, 



the flexoelectric contribution at the boundary surfaces is rep- 
resented by the last term of Eq. d4j) where 8^ > is the sur- 
face flexoelectric constant. This contribution depends linearly 
on the gradients of Q, and becomes irrelevant in the strong 
anchoring regime (because the integral of the associated free 
energy density can be rewritten as a surface integral over the 
boundary). 

It is often convenient to express the free energy in terms of 
dimensionless quantities. This decreases the number of pa- 
rameters necessary to describe the system, providing a min- 
imal set on which phase behaviour and electric field effects 
depend. The parameters are 
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The reduced temperature x multiplies the quadratic term of 
the dimensionless bulk free energy, whereas the gradient free 
energy term is proportional to the chirality K. This parameter 
can be considered as a measure of the twist present in the sys- 
tem; it is indeed defined as a ratio between the bulk and the 
gradient free energies. Lastly, parameters Th l and *£% 19 mea- 



sure the strength of bulk and surface flexoelectric couplings. 
Also notice that in strong anchoring regime W is set to infinity 



and £/ 



0. 



The equation of motion for Q i: 
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(3, + 3-V)Q-S(W,Q)=rH 



(9) 



where T is a collective rotational diffusion constant. The first 
term of the left hand side of Eq. © is the material deriva- 
tive describing the usual time dependence from advection by 
a fluid with velocity u. For the tensor field Q, this is general- 
ized by the term— 



S(W,Q) = (^D + co)(Q + I/3) + (Q- 
2^(Q + I/3)Tr(QW) 
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where Tr denotes the tensorial trace, while D = (W + W r ) /2 
and CO = (W — W r )/2 are the symmetric and the anti- 
symmetric part of the velocity gradient tensor W a $ = d^u a - 
The presence of S(W,Q) is due to the fact that the order pa- 
rameter distribution can be both rotated and stretched by the 
fluid. This is a consequence of the rod-like geometry of the LC 
molecules. The constant Z, depends on the aspect ratio of the 



molecules forming a given liquid crystal and controls whether 
the director field is flow-aligning in shear flow > 0.6 for 
y = 3), creating a stable response, or flow-tumbling, which 
gives an unsteady response < 0.6 for y = 3). The term on 
the right-hand side of Eq. © describes the relaxation of the 
order parameter towards the minimum of the free energy. The 
molecular field H which provides the driving motion is given 
by 

H = -^ + (I/3)Tr^. (11) 

with I being the unit matrix and J the free energy. 

The three-dimensional fluid velocity, u, obeys the continu- 
ity equation, 

3,p + 3 a (p^ a )=0 (12) 
and the Navier-Stokes equation, 

p(3, + updp)u a = 9p(n a p) +Ti3p(3 a wp + 3pM a ), (13) 

where p is the fluid density and T| is an isotropic viscosity. 
Note that the form of this equation is similar to that for a sim- 
ple fluid. However, the details of the stress tensor, n a p, re- 
flects the additional complication of cholesteric liquid crystals 
hydrodynamics. This tensor is explicitly given by: 
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Here Po is an isotropic background pressure. 

2.2 Numerical set up and mapping to physical units 

In this section we briefly describe the numerical aspects of our 
work. 

The size of the simulation box is equal to L x = 32, = 
32, L z = 32 lattice sites. Note that the size of the unit cell 
which minimises the free energy is in general different from 
(and typically larger than) the cholesteric pitch. Because we 
can only simulate an integer number of periodic unit cells 
within a periodic simulation domain, in order to account for 
this extra degrees of freedom we employ the following pro- 
cedure. Keeping fixed the size of the periodic domain, we 
rescale the size of the simulation box by a factor called "red- 
shift' In simulations with periodic boundaries in all 
directions, we have followed the recipe proposed in Ref. — , 

f This name refers to the shift, towards longer wavelengths, of the BPI Bragg 
back reflection with respect to that of the cholesteric phase—-. In other words, 
the unit cell size of the blue phase is typically larger than the cholesteric half 
pitch. 



in which the redshift is calculated at every time step by mini- 
mizing the free energy. Since the free energy functional com- 
prises terms up to quadratic order in gradients of Q, a rescal- 
ing of the unit cell dimension, K —> K/r (with r the redshift), 
modifies the gradient terms, that are in turn rescaled by a fac- 
tor of r per derivative. The free energy can be rewritten as 
/(Q) = r 2 a(dQ 2 ) + rb(dQ) + c, and the rescaling factor r* 
that minimizes the free energy is r* = —b/2a. In particular 
we assume a single redshift governing all the three dimen- 
sions^. We found that the presence of redshift does not affect 
the switching dynamics of BPs cellsi^. In simulations with 
homeotropic or homogenous anchoring, the redshift was set 
to one to avoid any change in the size of the confining cell. 
According to the experience gained in the bulk simulations, 
allowing for a finite redshift can lead to quantitative, but not 
qualitative, differences in the dynamics. 

The equilibrium structures of BPI and BPII in the absence 
of any field were obtained by relaxing the free energy to its 
minimum value by numerically solving Eq. ©. Initial config- 
urations were set according to the approximate solutions given 
in Ref. 16 , where it was shown that the topological character of 
the equilibrium defect structure at high chirality (where the 
bulk and the external field free energy terms become negli- 
gible compared to the gradient free energy and an analyti- 
cal solution is possible) is retained by reducing the chirality. 
These solutions are taken as initial conditions for the dynamic 
equations 22 ©, (IT2l) and (TT31) that are numerically integrated 
by using an hybrid scheme in which standard lattice Boltz- 
mann methods 11 are used to update the Navier-Stokes equa- 
tion while Eq. © is solved by using finite difference scheme 
based on a predictor corrector algorithm. This approach has 
been already successfully tested for other systems such as bi- 
nary fluids^ - —, active nematic s 11 ! 28 or blue phases 12 . A typ- 
ical simulation reported in this work was run in parallel on 8 
nodes with MPI architecture, and required about 5 days to be 
completed. 

In our calculations we fixed r = 0.3, £ = 0.7 and T| = 1.333 
as in previous numerical worksi2i2£. We have also set Ao = 
0.00308, K = 0.005, y = 3.775 for BPI, and A = 0.00075, 
K = 0.005, y = 3.086 for BPII. This choice guarantees we are 
in the correct region of the phase diagram 22 . In order to get 
from simulations to physical units we follow the approach of 
Refs.^2a2!i. Our parameters correspond to a blue phase with 
length scale fixed by the cholesteric pitch po, typically in the 
range of 100 — 500nm 16 , and (Frank) elastic constants equal 
to ~ 10 pN and ~ 16 pN for BPI and BPII respectively, as for 
a typical liquid crystalline material. One space and time LB 
unit correspond to 0.0125 ^m and 0.0013 jus for BPI, and to 
0.0125 ^m and 0.0017 jus for BPII, with a rotational viscos- 
ity equal to 1 Poise. The mapping of flexoelectric couplings 
can be done by matching the previously defined dimension- 
less quantities, and *E s fl , when defined in simulation and 



in physical units. For instance, for BPII a value of £y Z =0.3 
(which is in the middle of the range used in our runs) cor- 
responds to a system with (Frank) elastic constants equal to 
~ 10 pN as before, with a value of ~ 10 pC/m, subject to a 
field of ~ 15 Nljim. These values are typical of liquid crystal 
devices; to directly apply our results to this situation, we fur- 
ther need to select a liquid crystal for which dielectric effects 
are negligible, which is challenging at these field strengths - 
a possibility is to use bimesogenic materials which have near- 
zero dielectric anisotropy. On the other hand, a value of the di- 
mensionless surface flexoelectric coupling of 2^ =0.15 may 
correspond to the same BPII material, when the anchoring at 
the surface is W ~ 0.0007 N/m. For BPI, a similar mapping 
applies. 

Finally, to locate the defects in the simulation domain, we 
used the approach described in 12 , where a lattice point is con- 
sidered to belong to a defect if the local order at that point 
assumes values below a given threshold. Here the threshold is 
fixed to be 60% of the largest eigenvalue of Q in the bulk (i.e. 
away from all defects). 

3 Results 

3.1 Flexoelectric effects in blue phase samples with peri- 
odic boundary conditions 

We now discuss the results obtained by looking at the switch- 
ing dynamics of BPI and BPII unit cells with periodic bound- 
ary conditions. We will focus on the deformations of the de- 
fect structure observed when a field is applied to the sample 
along the z-direction ([0,0, 1] direction), only considering the 
effect of bulk flexoelectric couplings. The surface flexoelec- 
tric contribution vanishes with these boundary conditions, and 
the dielectric contribution is disregarded to focus on flexolec- 
tricity. For comparison, a single case in which the field is 
along the diagonal of the unit cell ([1,1,1] direction) will be 
also considered, to get a feel for the anisotropy of the response. 

BPs are arranged in a L x = L y = L z = 32 periodic unit cell. 
The initial equilibrium defect structures (see Fig{T]2 for BPI 
and Fig for BPII) are obtained by minimizing the free en- 
ergy by using exclusively the evolution equation Eq. © i2a2^. 
Notice that in this case the minimization is driven exclusively 
by the molecular field H. Unless otherwise stated, hydrody- 
namic interactions are switched on after this equilibration, and 
remain in place throughout the dynamical switching simula- 
tions. 

FigQ] shows the response of a periodic BPI unit cell un- 
der a bulk electric field applied along the z-direction. The 
field strength chosen c± 0.5) is large enough to modify 
the initial defect structure in an appreciable way. When the 
field is switched on, disclinations start to bend and twist, al- 
though the typical BPI defect structure is still recognizable 



(Fig 13?). As time goes on these defect distortions become 
more and more pronounced resulting into strongly bent arcs 
that span the entire unit cell and are aligned along the direc- 
tion of the field (FigHJ:). Later arcs reorganize to form "X- 
like"defect junctions (FigHJ/) that successively unbind at their 
center (FigQ]?) giving rise to columnar-like disclinations par- 
allel to the z-direction (FigQ]/*). FigQ] shows that the flexo- 
electric switching dynamics under a strong field is extremely 
different from previously studied cases in which only dielec- 
tric effects have been consideredJ^S. In the latter case the 
disclinations twist (for small field) or reconnect to lead to an 
amorphous state, in stark contrast with the regular structures 
of parallel defect lines obtained here. For lower fields (pro- 
vided these are still high enough to break the defect structure), 
the field-induced states are characterized by the same defect 
network, with the typical hexagonal pattern for the director 
profile. Interestingly, for those cases, during the switching on 
dynamics, double helical disclinations are observed initially, 
perpendicularly to the direction of the applied field - these 
later on merge and reorganize into the column-like pattern. 
Similar textures have already been described in Ref . 23 as in- 
termediate field-induced states, while in our simulations these 
are only transient defect structures. These differences might 
be due to different values of bulk and elastic constants consid- 
ered (for example, the theoretical description in 23 relaxes the 
one elastic constant approximation). 

A crucial point in the physics of devices is the understand- 
ing of the switching off dynamics, i.e. the time evolution of the 
system when the field is switched off. In particular one would 
like to answer the following question: Is the device switchable 
(i.e. does it recover the initial defect structure when the field 
is removed), or does it get stuck into a metastable state char- 
acterized by a different defect network ? In our experiment, as 
the flexoelectric coupling is switched off, the disclination lines 
start to bend, approach one another (Fig{T^), and merge into a 
connected structure that is similar in shape to the "X-like" state 
observed during the switching on dynamics but with thinner 
disclinations (FigUJi). Later on defects reorganize themselves 
into arcs that first span the entire system (Fig{TJ), then grad- 
ually stretch (FigQ]/) and finally split forming long undulated 
defect lines (Fig 03: and /), some of them in proximity of the 
center of unit cell, others closer to the boundaries. 

As Fig. 1 clearly shows, the system cannot recover the ini- 
tial state of FigQ]?, but it gets stuck into a metastable state (see 
Fig- 03) whose free energy is higher than the initial one. On 
the other hand this design can be still considered as a switch- 
able device since further simulations show that there exists a 
reversible cycle between the field-induced state in Fig. If and 
the metastable state in Fig. 1/. It is interesting to look also 
at the director profile of the state of Figfl]f induced by the 
flexoelectric coupling. In Fig 12 this is shown for a cross sec- 
tion in the xy plane taken at L z /2\ the structure consists of an 



hexagonal lattice of defects of — 1 / 2 topological charge that 
separate regions of splay-bend distortions around a line pass- 
ing through the center of the hexagon (i.e., the disclination 
in three dimensions). This spatial organization of the defects 
is in agreement with what found in Refs. 19 i 31 , where a similar 
structure has been identified as a stable configuration of the di- 
rector field, for either a spontaneously splayed nematic liquid 
crystal or a blue phase under an intense electric field, respec- 
tively. We find here that this state is naturally found in blue 
phases with flexoelectricity, and hence should be particularly 
relevant for the stabilised BPs recently found experimentally 
in Ref A 

A different kinetic pathway is observed when a field (still 
along the z-direction with £^ ~ 0.5) is applied on a BPII 
structure with periodic boundary conditions (see Figj3}. With 
respect to the BPI case the effect of the field is less dramatic 
and the typical (zero field) BPII network of defects is still 
recognizable even when the field is on. The persistence of 
a structure resembling (or topologically connected with) that 
of zero field, arguably suggests that the device should switch 
back to BPII when the field is taken off. This is indeed what 
happens. More precisely, when the field is switched on, the 
defect structure is initially bent and twisted at its center, while 
the regions in which the director field is distorted around a de- 
fect line widens (Fig 13?). The latter effect is made apparent by 
looking at the effective thickness of the disclination tubes. The 
device is driven towards a field induced state characterised by 
a defect structure similar to the zero-field (equilibrium) one 
(the shift along the z-direction is immaterial in view of the 
periodic boundary conditions, see FigOt). When the field is 
switched off, the disclination profile remains essentially undis- 
torted, with only a progressive slow reduction of the thickness 
of the defects (FigQ}/ and e) until the typical zero-field BPII 
structure (shifted along the z-direction) is recovered. Notice 
that only one cycle is here necessary to restore the initial de- 
fect pattern. 

We now look at the dependence of the defect dynamics with 
respect to the direction of the applied field. In FigH]we re- 
port the time evolution of the defect network of BPI when a 
field is applied along the diagonal of the unit cell direction 
([1,1,1]) which is not related to [0, 0, 1] by the cubic symmetry 
of the structure. Similarly to the case with the field along the 
z-direction, the strength of the applied field is large enough to 
significantly rewire the defect network. The steady-state field- 
induced state still has the defect lines broadly aligned along 
the field, but not as markedly as in the [0, 0, 1] case. When the 
field is switched off, the device gets again stuck in a metastable 
configuration, different than the one obtained when removing 
the field along the [0,0, 1] direction. 

Intriguingly, we find the switching-on dynamics depends 
strongly on field direction, and in all cases proceeds via multi- 
ple steps visible in the transient plateaus observed in the plots 




Fig. 1 Defect dynamics of the BPI phase during the switching on (a) — /)) -off (g) — /)) in a simulation box of size 32 x 32 x 32. Periodic 
boundary conditions are set in all directions. The electric field, applied along the z-direction ([0,0, 1]), is £^ ~ 0.5. It is switched off at 

t = 100 x 10 4 . One cycle is not enough to have a switchable device. A further application of the field reveals the device is switchable between 
the states /) and /). 



Fig. 2 Director profile on the xy plane taken at L z /2 of the metastable state /) of FigU Red circles indicate the presence of a charge —1/2 
topological defect. 




Fig. 3 Switching defect dynamics of a BPII phase of unit cell size 32x32x32 with periodic boundary conditions along all directions. The 
electric field, applied along the z-direction ([0,0, 1]), is ~ 0.5 (from a) to c)). It is switched off at t = 80 x 10 4 (from d) to /)). After one 
cycle, the system recovers the initial disclination structure (shifted along the z-direction), making the system switchable. 




Fig. 4 Evolution of the disclination network of the BPI phase during the switching on (a) —c)) -off (d) — /)) dynamics in a simulation box of 
size 32 x 32 x 32 with periodic boundary conditions in all directions. The electric field, applied along the diagonal direction ([1,1, 1]), is 
£^ ~ 0.5. It is switched off at t = 60 x 10 5 . When the field is on, the device is driven towards a metastable state. When the field is removed 
again a metastable state feauturing twisted disclinations is attained. 
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Fig. 5 Time evolution of the free energy in a BPI unit cell during a switching on-off dynamics with an electric field applied along the 
z-direction ([0,0, 1]) (plot (a)) and along the diagonal direction ([1,1, 1]) (plot (b)). The field is < Ej l ~ 0.5 in both cases. Interestingly the time 
evolution of the diagonal field case is ~ 6 time longer than the case in which the field is along the z-direction. (These free energy curves 
correspond to runs without hydrodynamics. This omission allows us to reach a larger value of the electric field, which emphasises the 
differences in the dynamics, avoiding at the same time numerical stability issues). 



of the free energy versus time (see Fig. |5}. 
3.2 BP samples under confinement 

We now turn to the switching dynamics of a field in confined 
BP samples. These are particularly relevant to real- world de- 
vices, as these necessarily include boundaries, where the di- 
rector field can easily be anchored, either chemically or me- 
chanically via rubbing - and the anchoring strength can often 
be controlled. In this Section we will only focus on BPI cells, 
which, as shown in the previous Section, lead to a more inter- 
esting dynamics in the bulk. We consider either homeotropic 
(normal) or homogeneous (planar) anchoring on upper and 
lower walls, and periodic boundary conditions in other direc- 
tions. To better compare with the results discussed previously, 
the size of the simulation box is set to 32 x 32 x 32. The field 
is then applied along the z-direction ([0,0,1]). Once more, 
its magnitude is strong enough that the initial defect network 
breaks up. The case of a surface flexoelectric coupling is also 
explored below - in this case the field effects are confined to 
the boundaries. 

The presence of confining surfaces strongly affects the sta- 
bility of defects in liquid-crystalline systems. In previous 
numerical studie s 32 i 33 , the existence of several defect struc- 
tures (sometimes described as "exotic" structures), very dif- 
ferent from that of the equilibrium BPI cell, have been pro- 
posed. A parallel array of double-helical disclinations, and 
two sets of undulating defect lines, were observed in a cell 
with homeotropic anchoring on both (parallel and flat) walls, 



when their distance was of the order of the helical pitch—. On 
the other hand, ring defects structures emerge when homoge- 
neous anchoring is set on both walls^. A key question, cru- 
cial for devices, concerns how to dynamically switch among 
different metastable states to kinetically reach, in a controlled 
way, a target structure. Recent studies 12 on the switching de- 
fect dynamics showed that applied electric fields could repre- 
sent an interesting route to reveal the existence of even more 
metastable states. A simple switching on-off schedule can 
lead, for example, to a bistable defect network 14 , where each 
state is uniquely selected by the application of the electric field 
along an appropriate direction. Interestingly, some key ques- 
tions arise when flexoelectric in confined samples is taken into 
account. What is the switching dynamics? Is this affected by 
the choice of different boundary conditions? 

To address these issues, in Figj6]we show the flexoelectric 
switching of a BPI cell with homeotropic anchoring (with a 
slight pretilt of ~ 10° from the z-direction) on both walls. We 
have also set £y Z ~ 1 (applied along the z-direction) and Wo = 
0.1, which still corresponds to moderately strong anchoring. 

The first thing to note is that the equilibrium defect structure 
(FigJSJO is different from the typical one observed for the BPI 
at equilibrium with periodic boundaries, although the basic 
BPI topology is still recognisable (for instance in the discli- 
nations in the interior of the cell). This metastable structure 
is the same as that reported in Refs.J^*^. A second, key, ob- 
servation is that the dynamics of flexoelectric switching are 
significantly affected by the anchoring. Thus, the normal an- 



choring conditions in Fig.[6]lead to a visible initial bend of the 
diagonal disclination. This is followed later on by the forma- 
tion of straight defect lines as in the bulk case (see Figj6]b and 
c), which branch out at the boundary to end up parallel to the 
confining planes there (see Fig.[6]i-/). Once more, upon field 
removal the system drifts away from the field-induced state in 
Fig. [6f into another metastable state, featuring undulating de- 
fect lines close to but distinct from those observed in the bulk 
(compare Fig. [6} with Fig.[TJ). 

Homogeneous boundary conditions, as can be achieved by 
rubbing the liquid crystal at the boundaries, lead to yet differ- 
ent structures in equilibrium, and to new kinetic pathways dur- 
ing switching on and off a field. The behaviour of a BPI cell 
with such planar anchoring (along the x direction) is shown in 
Fig. [71 The flexoelectric coupling is again given by Tj t ~ 1, 
and the anchoring Wo = 0.1 as in Fig. [6] The equilibrium de- 
fect network is now made up of an array of isolated twisted 
ring disclinations located in the bulk of the device (Fig|7j2). 
These structures have already been discussed in Ref . 33 and, 
as a field induced state, in the bistable BP device proposed in 
Ref.—: we do not dwell on them here, and rather focus on the 
field-induced aspects. When the field is switched on, the net- 
work is initially bent and stretched (Fig|3? and c) with discli- 
nations again alinging parallel to the field until they touch the 
walls. Later on they recombine into elongated rings now con- 
nected to defects pinned up at the walls ((FigJTJ/, e and f). 
As in the previous cases, when the field is removed defects 
depin from the surfaces and form connected ring arcs in the 
bulk (from Fig|7^ to j). Then the network develops branching 
points which are unstable and annihilate in the bulk leaving a 
state with arcs of defects spanning the whole device (Fig|7](: 
and /). 

Lastly, we investigate the switching dynamics in response to 
surface flexoelectricity. For strong anchoring these contribu- 
tions can be rewritten as surface terms, so one may expect that 
their presence may only lead to deviations from the equilib- 
rium bulk disclinations in strongly confined samples, such as 
the ones we are considering. As we mentioned previously, be- 
cause the field couples with the derivative of the order param- 
eter (see Eq|4]), in the case of strong anchoring the coupling is 
identically zero at the boundaries and the surface flexolectric 
term is not effective. However for finite anchoring (as the case 
of Wo = 0.1 we consider), this is not true anymore, and the 
director can deviate from the anchoring direction, at the cost 
of some surface energy. 

To address the role of the surface flexoelectric coupling, we 
simulated the evolution of disclinations in a confined BPI cell 
with homogeneous anchoring (FigHJ). The magnitude of the 
field, applied along the z-direction, is such that 2^ ~ 0.15 
(a pretilt of ~ 10° from the plane is also included). Now 
the switching leads to a shrinkage and eventual annihilitation 
of the twisted disclination loops (the full kinetic pathway is 



shown in FigH^-i). The final state is defect- free (Fig.©, and 
the associated director profile has a lot of splay-bend defor- 
mations. Intriguingly, this state is similar to the spontaneously 
splayed nematics in an achiral liquid crystal with a strong flex- 
oelectric coupling, predicted in—. As expected, the state re- 
mains disclination free even when the field is removed. 

4 Conclusions 

To summarise, we have numerically investigated the switching 
dynamics of cholesteric blue phases in the presence of a flex- 
oelectric coupling with an electric field. The choice of bound- 
ary conditions (periodic boundaries, or confined samples with 
homogeneous or homeotropic anchoring on both walls) and 
of the dynamical schedule of the field application crucially 
affect the dynamics, yielding several intriguing defect struc- 
tures. Some of these resemble the zero-field equilibrium net- 
work (e.g. for BPII), while others organize forming novel 
field induced stable or metastable states. The most frequently 
observed structure is one where the disclinations are aligned 
(more or less markedly, according to the different cases) along 
the field direction. 

With periodic boundary conditions (corresponding to the 
bulk of a large sample) and with a field applied along the 
[0,0, 1] direction, BPI defects align along the direction of the 
applied field via a complex kinetic pathway which includes 
the formation of intermediate "X-like"disclination junctions. 
Once the field is switched off, the network is stuck in a zero- 
field metastable state. Remarkably, a further application of 
the field restores the field induced states, thereby making 
these cells switchable. When the field is applied along the 
[1,1,1] direction, the alignment is imperfect and the dynamics 
is slower. 

Importantly, in both cases, as for most of the other geome- 
try/boundary condition choices, upon switching off, the field- 
induced BPI structure drifts away from the field-induced con- 
formation and gets stuck into a metastable phase, the details 
of which depend on boundary conditions, field directions etc.. 
This means that, while the field-aligned disclination state is 
switchable, it cannot, at least under the conditions explored by 
our simulations, be used as the basis for multistable energy - 
saving devices such as those described in—. (However we 
have laid the formulation for future computational searches 
for these). The inability to reform an equilibrium BPI struc- 
ture also suggests that large memory effects are present. 

Confined samples lead to qualitatively different results. 
With normal, or homeotropic, anchoring at the boundaries, 
and bulk flexoelectric coupling only, the field-aligned discli- 
nations branch at or close to the wall to end up parallel to the 
boundaries there. With homogeneous anchoring the defects 
without a field form twisted rings; these elongate and align 
with the field when bulk flexoelectricity is switched on. Sur- 




Fig. 6 Evolution of the BPI disclination network in a confined sample of size 32 x 32 x 32 with homeotropic anchoring. Top and bottom walls 
(in grey) are oriented in the xy plane. The switching on is from a) to /) and the switching off is from g) to i). The field, applied along the 
z-direction ([0,0, 1]), is such that £^ ~ 0.5. It is switched off at t = 100 x 10 4 . The switching on drives the device towards a field induced 
state /) in which column-like defects are connected at the walls, while, when the field is removed, the cell gets stuck into a metastable state 
with straight disclinations close to the walls coexisting with undulating lines in the middle of the BPI cell. 




Fig. 7 Evolution of the BPI disclination network in a confined sample of size 32 x 32 x 32 with homogeneous anchoring. Top and bottom 
walls (in grey) are oriented in the xy plane. The switching on is from a) to /) and the switching off is from g) to /). The field, applied along 
the z-direction ([0,0, 1]), is such that £^ ~ 1. It is switched off at t = 128 x 10 4 . Twisted ring defects are recognisable during the switching 
on, while a novel metastable state /), with arcs of defects spanning the device, is obtained when the field is removed. 




Fig. 8 Evolution of the BPI disclination network in a confined sample of size 32 x 32 x 32 with homogeneous anchoring and with surface 
flexoelectric coupling. Top and bottom walls (in grey) are along xy plane. The switching on is from a) to i), where all defects disappear. The 
field, applied along the z-direction ([0,0, 1]), is such that *E% ~ 0.15. 



X 



Fig. 9 3D director profile of the field-induced defect-free state in Fig.[8j. 



face flexoelectric contributions, relevant at low to intermediate 
anchoring strength, can drive the device into yet more configu- 
rations, such as a strongly splay ed-bent nematic conformation 
(see Fig. [9]). 

Our results give, we believe, a convincing demonstration 
that flexoelectricity and geometry are important factors in se- 
lecting defect networks in blue phases. Given the importance 
of flexoelectric coefficients in the materials used to build the 
new blue phases used in Ref.-, we hope our results will be 
useful when designing new blue phase-based liquid crystal de- 
vices. 



Acknowledgements 

We thank EPSRC EP/J007404 for funding. MEC holds a 
Royal Society Research Professorship. 



References 

1 H. Kikuchi, M. Yokota, Y. Hisakado, H. Yang and T. Ka- 
jiyama, Nat. Mater. 1 64 (2002). 



2 Y. Hisakado, H. Kikuchi, T. Nagamura and T. Kajiyama, 
Adv. Mat. 17 96 (2004). 

3 H. J. Coles and M. N. Pivnenko, Nature 436, 997 (2005). 

4 F. Castles, F. V. Day, S. M. Morris, D-H. Ko, D. J. Gar- 
diner, M. M. Qasim, S. Nosheen, P. J. W. Hands, S. S. 
Choi, R. H. Friend and H. J. Coles, Nat. Mat. 11, 599 
(2012). 

5 E. Karatairi, B. Rozic, Z. Kutnjak, V. Tzitzis, G. Nounesis, 
Phys. Rev. E, 81 041703 (2010). 

6 www.physorg.com/newsl29997960.html (2008). 

7 C. Denniston and J. M. Yeomans, Phys. Rev. Lett. 87, 
275505 (2001); A. Tiribocchi, G. Gonnella, D. Maren- 
duzzo, and E. Orlandini Appl. Phys. Lett. 97, 143505 
(2010). 

8 H. J. Coles, M. J. Clarke, S. M. Morris, B. J. Broughton, 
and A. E. Blatch, /. Appl. Phys. 99, 34104 (2006). 

9 F. Castles, S. M. Morris, E. M. Terentjev, and H. J. Coles, 
Phys. Rev. Lett. 104, 157801 (2010). 

10 O. Henrich, D. Marenduzzo, K. Stratford, M. E. Cates, 
Phys. Rev. E 81, 031706 (2010). 

11 M. E. Cates, O. Henrich, D. Marenduzzo, K. Stratford, 
Soft Matter 5, 3791-3800 (2009). 

12 A. Tiribocchi, G.Gonnella, D. Marenduzzo and E. Orlan- 



dini, Soft Matter 7, 3295 (2011). 

13 O. Henrich, K. Stratford, D. Marenduzzo, M. E. Cates, 
Proc. Natl Acad. Set USA 107, 13212 (2010); O. Henrich, 
K. Stratford, M. E. Cates, and D. Marenduzzo, Phys. Rev. 
Lett. 106, 107801 (2011). 

14 A. Tiribocchi, G.Gonnella, D. Marenduzzo, E. Orlandini, 
F. Salvadore, Phys. Rev. Lett. 107, 237803 (2011). 

15 A. N. Beris, B, J. Edwards, Thermodynamics of Flowing 
Systems, Oxford University Press, Oxford (1994). 

16 D. C. Wright, N. D. Mermin, Rev. Mod. Phys. 61, 385 
(1989) and references therein. 

17 P. G. de Gennes and J. Prost, The Physics of Liquid Crys- 
tals, 2nd Ed., Clarendon Press, Oxford (1993). 

18 S. Chandrasekhar, Liquid Crystals, Cambridge University 
Press, (1980). 

19 G. P. Alexander, J. M. Yeomans, Phys. Rev. Lett. 99, 
067801 (2007). 

20 D. Marenduzzo, E. Orlandini, J. M. Yeomans, Phys. Rev. 
Lett. 92, 188301 (2004). 

21 H. Grebel, R. M. Hornreich, and S. Shtrikman, Phys Rev. 
A 28, 1114(1983). 

22 A. Dupuis, D. Marenduzzo, J. M. Yeomans, Phys. Rev. E 
71,011703 (2005). 

23 G. P. Alexander, J. M. Yeomans, Liq. Cryst. 36, 1215-1227 
(2009). 

24 A. Dupuis, D. Marenduzzo, E. Orlandini, J. M. Yeomans, 
Phys. Rev. Lett. 95, 097801 (2005). 

25 A. G. Xu, G. Gonnella and A. Lamura, Physica A 362, 42 
(2006). 

26 A. Tiribocchi, N. Stella, G. Gonnella, A. Lamura, Phys. 
Rev. E 80, 026701 (2009). 

27 G. Gonnella, A. Lamura, A. Piscitelli and A. Tiribocchi, 
Phys. Rev. £82, 046302 (2010). 

28 D. Marenduzzo, E. Orlandini, M. E. Cates, J. M. Yeomans 
Phys. Rev. E76, 031921 (2007). 

29 G. P. Alexander, J. M. Yeomans, Phys. Rev. E. 74, 061706 
(2006). 

30 J. Fukuda, M. Yoneya, and H. Yokoyama, Phys. Rev. E 80, 
031706 (2009). 

31 R. M. Hornreich, M. Kugler, S. Shtrikman, Phys. Rev. Lett. 
54, 2099 (1985). 

32 J. Fukuda, S. Zumer, Phys. Rev. Lett. 104, 017801 (2010). 

33 J. Fukuda, S. Zumer, Phys. Rev. Lett. 106, 097801 (2011). 



